# Replication script for 
# Mędrzycki, Piotr; Jarzyna, Ingeborga ; Tokarska-Guzik, Barbara; Obidziński, Artur; Sotek, Zofia; Pabjanek, Piotr; Pytlarczyk, Adam; Sachajdakiewicz, Izabela, 2017, "Replication Data for: Simple yet effective: historical proximity variables improve the species distribution models for invasive giant hogweed (Heracleum mantegazzianum s.l.) in Poland", doi:10.15139/S3/BEBL6N, UNC Dataverse, DRAFT VERSION [UNF:6:0jm30/bymUtZd8iWfkgB1g==] 

# Adapted from the Biomod2 Multi_species_computation.R script with unvaluable help of Damien Georges <damien.georges2@gmail.com> 

rm(list=ls())

# 1. loading species occurrences data

# 1. loading species occurrences data

library(biomod2)

# 
# 
# DataSpecies <- read.csv( system.file( 
#   
#   "external/species/mammals_table.csv", 
#   
#   package="biomod2"))
# 
# 
# 
# head(DataSpecies)







###################################################

### code chunk number 3: LoadEnv_1

###################################################

# 2. loading environmental data



# Environmental variables extracted from WorldGrids 3rd level

require(raster)

myExplPL <- stack(
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/DEMSRE3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/SLPSRT3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/cntgad3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/etmnts3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/evmmod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/evsmod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g01esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g01igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g02esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g02igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g03esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g04esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g04igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g05esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g06esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g07esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g09esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g10esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g10igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g11igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g12esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g12igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g13esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g14esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g15esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g16esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g17esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g18esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g19esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g20esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g21esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/g22esa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gabhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gacgem3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gachws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/galhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/ganhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/garhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gathws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gchhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gclhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gcmhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gcrhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/geaisg3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gflhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gfrhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gglhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/ggyhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/ghshws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gkshws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/glcesa3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/glcjrc3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/glphws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/glvhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/glwwwf3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/glxhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gnthws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gphhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gplhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gpthws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gpzhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/grghws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gschws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gsnhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gsthws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gumhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/gvrhws3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/iflgre3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/inmsre3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/inssre3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l01igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l02igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l03igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l04igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l05igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l06igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l07igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l08igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l09igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l10igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l11igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l12igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l13igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l14igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l15igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l16igb3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/lasmod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/lmbgsh3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/lmtgsh3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/ln1dms3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/ln2dms3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/lnmdms3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/opisre3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/smkisr3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/tdhmod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/tdlmod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/tdmmod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/tdsmod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/tnmmod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/tnsmod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/twisre3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/tx1mod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/tx1mod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/tx1mod3a.tif"),
  raster("/home/piotrmed/GIS/WorldGrids/GDALed/l3pobi3b.tif")
)

# Import of 5 regions boundaries shapefile----------
Study_area_SHP <- shapefile("5Regions_4326.shp")
plot(Mask)
plot(Study_area_SHP, add=T)
# Converting mask to raster and reclassification to binary grid------------
Study_area_MASK <- rasterize(Study_area_SHP, 
                             myExplPL, 
                             fun='last',
                             background=NA,
                             mask=FALSE,
                             update=FALSE, 
                             updateValue='all',
                             filename="", 
                             getCover=F, 
                             silent=FALSE)


plot(Study_area_MASK)

Study_area_MASKNA <- reclassify(Study_area_MASK, c(-Inf,0.25,NA, 0.25,10,1))
plot(Study_area_MASKNA)
describe(Study_area_MASKNA)
Study_area_MASKNA <- na.omit(Study_area_MASKNA)


# Additional Presence, count and proximity parameters



HM_hist <- shapefile("Sites_2007_5REG_EPSG4326a.shp")
proj4string(HM_hist) <- CRS("+init=epsg:4326")
proj4string(Study_area_SHP) <- CRS("+init=epsg:4326")
HM_hist_SA_NA <- HM_hist
# Selecting the date-depending datasets for the SA only -----
HM_1985 <- subset(HM_hist_SA_NA, AGE_CLASS=="1985")
HM_2002 <- subset(HM_hist_SA_NA, AGE_CLASS<= "2002")
HM_2007 <- subset(HM_hist_SA_NA, AGE_CLASS<= "2007")



plot(HM_1985)
plot(Study_area_SHP, add=T)

# Calculating quantitative variables----
require(raster)
# Sites presence before ----
HM_1985_Count <- rasterize(HM_1985, 
                           Study_area_MASK, 
                           field='Abundance',
                           fun='count',
                           background=NA,
                           mask=FALSE,
                           update=FALSE, 
                           updateValue='all',
                           filename="", 
                           getCover=F, 
                           silent=FALSE)
HM_1985_Count[is.na(HM_1985_Count)] <- 0
hist(HM_1985_Count)
density(HM_1985_Count)

# 2002
HM_2002_Count <- rasterize(HM_2002, 
                           Study_area_MASK, 
                           field='Abundance',
                           fun='count',
                           background=NA,
                           mask=FALSE,
                           update=FALSE, 
                           updateValue='all',
                           filename="", 
                           getCover=F, 
                           silent=FALSE)
HM_2002_Count[is.na(HM_2002_Count)] <- 0
hist(HM_2002_Count)
density(HM_2002_Count)
# 2007
HM_2007_Count <- rasterize(HM_2007, 
                           Study_area_MASK, 
                           field='Abundance',
                           fun='count',
                           background=NA,
                           mask=FALSE,
                           update=FALSE, 
                           updateValue='all',
                           filename="", 
                           getCover=F, 
                           silent=FALSE)
HM_2007_Count[is.na(HM_2007_Count)] <- 0
hist(HM_2007_Count)
density(HM_2007_Count)
# option, when there are data for historic individual nuumber
# # Total abundance before 1985----
# HM_1985_Abund <- rasterize(HM_1985, 
#                           Study_area_MASK, 
#                           field='Abundance',
#                           fun='sum',
#                           background=NA,
#                           mask=FALSE,
#                           update=FALSE, 
#                           updateValue='all',
#                           filename="", 
#                           getCover=F, 
#                           silent=FALSE)
# Site presence/absence --------
HM_1985_01 <- rasterize(HM_1985, 
                        Study_area_MASK, 
                        field='Abundance',
                        fun='last',
                        background=NA,
                        mask=F,
                        update=FALSE, 
                        updateValue='all',
                        filename="", 
                        getCover=F, 
                        silent=FALSE)


HM_1985_01 <- reclassify(HM_1985_01,
                         c(-Inf,0.25,NA, 0.25,+Inf,1))
HM_1985_01 <- mask(HM_1985_01,Study_area_MASKNA)
HM_1985_01[is.na(HM_1985_01)] <- 0
plot(Study_area_MASKNA)
plot(HM_1985_01, add=T)
density(HM_1985_01)   

# 2002
HM_2002_01 <- rasterize(HM_2002, 
                        Study_area_MASK, 
                        field='Abundance',
                        fun='last',
                        background=NA,
                        mask=F,
                        update=FALSE, 
                        updateValue='all',
                        filename="", 
                        getCover=F, 
                        silent=FALSE)


HM_2002_01 <- reclassify(HM_2002_01,
                         c(-Inf,0.25,NA, 0.25,+Inf,1))
HM_2002_01 <- mask(HM_2002_01,Study_area_MASKNA)
HM_2002_01[is.na(HM_2002_01)] <- 0
plot(Study_area_MASKNA)
plot(HM_2002_01, add=T)
density(HM_2002_01)   
# 2007
HM_2007_01 <- rasterize(HM_2007, 
                        Study_area_MASK, 
                        field='Abundance',
                        fun='last',
                        background=NA,
                        mask=F,
                        update=FALSE, 
                        updateValue='all',
                        filename="", 
                        getCover=F, 
                        silent=FALSE)


HM_2007_01 <- reclassify(HM_2007_01,
                         c(-Inf,0.25,NA, 0.25,+Inf,1))
HM_2007_01 <- mask(HM_2007_01,Study_area_MASKNA)
HM_2007_01[is.na(HM_2007_01)] <- 0
plot(Study_area_MASKNA)
plot(HM_2007_01, add=T)
density(HM_2007_01)   
# Proximity distance to the nearest point -------
# 1985 
HM_1985_Prox <- distanceFromPoints(Study_area_MASKNA, 
                                   HM_1985) 

HM_1985_Prox[is.na(HM_1985_Prox)] <- 0
plot(HM_1985_Prox)

# 2002
HM_2002_Prox <- distanceFromPoints(Study_area_MASKNA, 
                                   HM_2002) 
plot(HM_2002_Prox)
HM_2002_ProxMask <- crop(HM_2002_Prox,
                         myExplPL,
                         inverse=T, 
                         maskvalue=1) 
HM_2002_Prox[is.na(HM_2002_Prox)] <- 0
plot()
plot(HM_2002_ProxMask)
density(Study_area_MASKNA)

# 2007 
HM_2007_Prox <- distanceFromPoints(Study_area_MASKNA, 
                                   HM_2007) 
plot(HM_2007_Prox)
HM_2007_ProxMask <- crop(HM_2007_Prox,
                         myExplPL,
                         inverse=T, 
                         maskvalue=1) 
HM_2007_Prox[is.na(HM_2007_Prox)] <- 0
plot()
plot(HM_2007_ProxMask)
density(Study_area_MASKNA)

# Import of data from 2012 and 2013 yr for all country----------
# Data already crppped to the range of 5 regions


HM_2012 <- shapefile("Sites_20120229_5REG_EPSG4326a.shp")



# Change to relative age and substituting NAs with 0 -----------

# Presence / Absence
HMD.25.1 <- HM_1985_01



HMD.10.1 <- HM_2002_01


HMD.5.1 <- HM_2007_01


# Site count

HMD.25.C <- HM_1985_Count


HMD.10.C <- HM_2002_Count


HMD.5.C <- HM_2007_Count


# Site proximity

HMD.25.P <- HM_1985_Prox


HMD.10.P <- HM_2002_Prox


HMD.5.P <- HM_2007_Prox


HMD.25.10.P <- stack(HMD.25.P, HMD.10.P)


HMD.25.10.5.P <- stack(HMD.25.P,HMD.10.P,HMD.5.P)
HMD.25.10.5.P <- stack(HMD.25.C,HMD.10.C,HMD.5.C)
HMD.25.10.5.1 <- stack(HMD.25.1,HMD.10.1,HMD.5.1)



HMD.25.CP <- stack(HMD.25.P, HMD.25.C)
HMD.25.10.CP <- stack(HMD.25.P, HMD.25.C, HMD.10.P, HMD.10.C)

HMD.25.10.C <- stack(HMD.25.C,HMD.10.C)

HMD.25.10.5.CP <- stack(HMD.25.P, HMD.25.C, HMD.10.P, HMD.10.C, HMD.5.P, HMD.5.C)

HMD.5.CP <- stack(HMD.5.P, HMD.5.C)
HMD.10.CP <- stack(HMD.10.P, HMD.10.C)
HMD.25.10.5.P <- stack(HMD.25.P,HMD.10.P,HMD.5.P)
HMD.25.10.5.C <- stack(HMD.25.C,HMD.10.C,HMD.5.C)
HMD.25.10.5.1 <- stack(HMD.25.1,HMD.10.1,HMD.5.1)
HMD.25.10.C <- stack(HMD.25.C,HMD.10.C)

HMD.0 <- Study_area_MASK

HMD.10.5.1 <- stack(HMD.10.1, HMD.5.1)
HMD.10.5.C <- stack(HMD.10.C, HMD.5.C)
HMD.10.5.P <- stack(HMD.10.P, HMD.5.P)
HMD.10.5.CP <- stack(HMD.10.CP, HMD.5.CP)


# Checking distribution
hist(HM.5.C)
hist(HM.5.P)
hist(HM.5.1)

hist(HM.10.C)
hist(HM.10.P)
hist(HM.10.1)

hist(HM.25.C)
hist(HM.25.P)
hist(HM.25.1)


###################################################

### code chunk number 4: Loop_1

###################################################

# define the species of interest

# Setting the reposnse variables


sp.name <- "HM"

# myResp = as.numeric(DataSpecies[,sp.name])

myResp <- SpatialPoints(HM_2012)
plot(myResp, add=T)

myRespName = sp.name



# I define the list of datasets

ls(pattern="^HMD.")

#hab = list.files(getwd(), pattern="HM*$", full.names=FALSE) 

myExplCPNames <- ls(pattern="^HMD.")
myExplCPNamesAdd <- ls(pattern="^HMD.10.5.")



# It works!!
for(j in myExplCPNamesAdd){
  
  nam <- paste("myExplD",j, sep=".")
  assign(nam, stack(get(j),myExplPL))
}


myExplD.HMD.0 <- myExplPL
#rm(list=ls(pattern="myExplD.HM"))


myExplPOOL <- ls(pattern="myExplD.HMD.")
myExplPOOLAdd <- c("myExplD.HMD.5.CP","myExplD.HMD.10.CP" , "myExplD.HMD.25.10.5.P","myExplD.HMD.25.10.5.C","myExplD.HMD.25.10.5.1","myExplD.HMD.25.10.C")
myExplPOOLAdd <-
  
  # loop on datasets instead of species == applying the same functions to each
  
  
  
  for(myExpl.n in myExplPOOLAdd){
    
    # I want to name the models after the dataset chosen
    
    myExplName = myExpl.n
    
    
    
    cat('\n',myRespName ,"with", myExplName,'modeling...')  
    
    ### definition of data 
    
    ## i.e keep only the column of our species
    
    # here it is impossible to adapt
    
    myExpl <- get(myExpl.n)
    
    
    
    ### Initialisation
    library(biomod2) 
    myBiomodData <- BIOMOD_FormatingData(resp.var = myResp,
                                         expl.var = myExpl,
                                         resp.xy = NULL,
                                         resp.name = myRespName,
                                         PA.nb.rep = 5,
                                         PA.nb.absences = 1000,
                                         PA.strategy = 'random',
                                         PA.dist.min = 0,
                                         PA.dist.max = NULL,
                                         PA.sre.quant = 0.025,
                                         PA.table = NULL,
                                         na.rm = TRUE)
    
    
    
    
    
    ### Options definition
    
    myBiomodOption <- BIOMOD_ModelingOptions()
    
    
    
    
    ### Modelling 
    
    myBiomodModelOut <- BIOMOD_Modeling( 
      myBiomodData, 
      models = c('SRE','CTA','RF','MARS','FDA', 'GBM', 'ANN'), 
      models.options = myBiomodOption, 
      NbRunEval=3, 
      DataSplit=80, 
      Prevalence=0.5, 
      VarImport=3,
      models.eval.meth = c('TSS','ROC'),
      SaveObj = TRUE,
      rescal.all.models = TRUE,
      do.full.models = FALSE,
      modeling.id = paste("with_",myExplName,sep=""))
    
    
    
    ### save models evaluation scores and variables importance on hard drive
    
    modeling.id = paste("with_",myExplName,sep="")
    
    capture.output(get_evaluations(myBiomodModelOut),
                   file=file.path(myRespName, 
                                  paste(myExplName,"_formal_models_evaluation.txt", sep="")))
    
    
    capture.output(get_variables_importance(myBiomodModelOut),
                   file=file.path(myRespName,
                                  paste(myExplName,"_formal_models_variables_importance.txt",
                                        sep="")))
    
    
    ### Building ensemble-models
    
    myBiomodEM <- BIOMOD_EnsembleModeling( 
      modeling.output = myBiomodModelOut,
      chosen.models = 'all',
      em.by='all',
      eval.metric = c('TSS'),
      eval.metric.quality.threshold = c(0.7),
      prob.mean = T,
      prob.cv = T,
      prob.ci = T,
      prob.ci.alpha = 0.05,
      prob.median = T,
      committee.averaging = T,
      prob.mean.weight = T,
      prob.mean.weight.decay = 'proportional' )
    
    
    capture.output(get_evaluations(myBiomodEM),
                   file=file.path(myRespName, 
                                  paste(myExplName,"_ensemble_models_evaluation.txt", sep="")))
    
    
    
    
    capture.output(get_variables_importance(myBiomodEM),
                   
                   file=file.path(myRespName, 
                                  
                                  
                                  
                                  
                                  paste(myExplName,"_ensemble_models_variables_importance.txt",
                                        sep="")))
    
    
    
    
    ### Make projections on current variable
    
    myBiomodProj <- BIOMOD_Projection(
      
      modeling.output = myBiomodModelOut,
      
      new.env = myExpl,
      
      proj.name = 'current',
      
      selected.models = 'all',
      
      binary.meth = 'TSS',
      
      compress = 'xz',
      
      clamping.mask = F,
      
      output.format = '.grd')
    
    
    ### Make ensemble-models projections on current variable
    
    myBiomodEF <- BIOMOD_EnsembleForecasting( 
      
      projection.output = myBiomodProj,
      
      EM.output = myBiomodEM)
    
  }


# Metamodelling of model results------------

# Non standardized model 

Data_for_metamodels <- read.csv("Data_for_metamodels.csv", sep=";")
library(Hmisc)
describe(Data_for_metamodels[c(4,7,10,13:16)])
library(latticist)
latticist(Data_for_metamodels)
names(Data_for_metamodels[c(4,7,10,13:16)])
library(Boruta)
TSS_bor <- Boruta(TSS~., data=Data_for_metamodels[c(4,7,10,13:16)])
attStats(TSS_bor)
plot(TSS_bor)

# TSS standardized on model type
## TSS standardized inside algorithm
## Model with TSS only

TSS_norm_bor <- Boruta(TSS_norm~., data=Data_for_metamodels[c(4,7,10,13:16)])
attStats(TSS_norm_bor)
plot(TSS_norm_bor)

# Response plots with looping (for additional visualisation) ---------------------------------------


# Looping over variables ranked by importance 
set.seed(135)
library(randomForest)
TSS_rf <- randomForest(TSS~., data=Data_for_metamodels[c(4,7,10,13:16)])
TSS_imp <- importance(TSS_rf)
TSS_impvar <- rownames(TSS_imp)[order(TSS_imp[, 1], decreasing=TRUE)]

op <- par(mfrow=c(2, 3), las=2)

op2 <- par(mfrow=c(1, 1), las=2)
for (i in seq_along(TSS_impvar)) {
  partialPlot(TSS_rf, Data_for_metamodels[c(4,7,10,13:16)], 
              TSS_impvar[i], xlab=TSS_impvar[i],
              main=paste("TSS Raw ~", TSS_impvar[i]),
              ylim=c(0.01,40), cex.axis=0.6, cex.lab=0.7, cex.main=0.7)  	 
}


